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Abstract 

Zirconium alloys are used as nuclear fuel cladding material due to their mechanical 
and corrosion resistant properties together with their favorable cross-section for neutron 
scattering. At running conditions, however, there will be an increase of hydrogen in 
the vicinity of the cladding surface at the water side of the fuel. The hydrogen will 
diffuse into the cladding material and at certain conditions, such as lower temperatures 
and external load, hydrides will precipitate out in the material and cause well known 
embrittlement, blistering and other unwanted effects. Using phase-field methods it is 
now possible to model precipitation build-up in metals, for example as a function of 
hydrogen concentration, temperature and external load, but the technique relies on input 
of parameters, such as the formation energy of the hydrides and matrix. To that end, we 
have computed, using the density functional theory (DFT) code GPAW, the latent heat 
of fusion as well as solved the crystal structure for three zirconium hydride polymorphs: 
<5-ZrHi.6, 7-ZrH, and e-ZrH2. 

Introduction 

The presence of hydrogen in certain metals such as Hf-, Ti- and Zr-based alloys can be 
a major cause of embrittlement due to precipitation of hydrides (MH^) pQ. There are at 
least three mechanisms associated with the hydrogen- induced mechanical degradation 
of metals [2]. These are (1) phase transformations, such as hydride precipitation induced 
by combined presence of H and stress, (2) H-enhanced local plasticity, and (3) weakening 
of grain boundaries by H. In the following we will take a closer look at the first of these. 

Hydrogen embrittlement may occur by repeated hydride precipitation and cleavage 
on front of a growing crack. The phenomena involve hydrogen transport and phase trans- 
formation (nucleation and growth) near a crack tip which is affected by the combined 
effects of stress concentration, temperature distribution and accommodation modes be- 
tween hydride and matrix. The situation in Zr-base systems is more complex by the 
existence of different phases; the hexagonal close-packed (hep) low temperature a-phase, 
and the body-centered cubic (bec) high temperature /3-phase, which have different sus- 
ceptibilities to hydrogen embrittlement. The a-phase being more susceptible and having 
a lower hydrogen solubility and diffusivity than the /3-phase. Alloys containing a mix- 
ture of a- and /3-phase usually crack at the a//3-phase interface J3J. In presence of flaws, 
notches or other stress concentrators, the stress-directed movement of hydrogen may 



also result in a sub-critical crack growth mechanism, known as delayed hydride cracking 
(DHC). Under DHC, hydrogen in solid solution diffuses to the high-stress region in front 
of the flaw, where it precipitates as brittle hydrides. The hydrided region grows larger 
and denser as more hydrogen arrives to the tip of the flaw, until the material is no longer 
able to withstand the local tensile stress. Brittle or semi-brittle fracture then ensues, 
and the flaw grows by an increment that is comparable in size to the hydrided region 
ahead of the flaw tip. The process then repeats itself, thus leading to stepwise crack 
growth that leaves characteristic striations on the fracture surfaces [I]. 

It is today possible to model general hydride micro-structure evolution in a primary 
matrix phase using a phase-field model[5j. For Zr in particular, Ma et al. used a phase- 
field kinetic model in a two-dimensional setting to investigate the hydride morphology 
as a function of external load 0, Ej. In any case, the usage of phase-field methods needs 
input of parameters, such as the latent heat of formation, and the accuracy of the model 
is determined by the source and precision of these parameters. However, only a few of 
the needed parameters that have been published are from reliable experimental sources 
and some ambiguities still exist regarding, for example, the atomic structure, elastic 
coefficients, and phase-diagram of the H-Zr system. The atomic structure used today 
are from neutron experiments that dates back some 50 years and it is quite debatable 
whether one should trust them fully. A number of recent works have used ab-initio 
computational methods to try to improve the picture somewhat and several groups have 
used DFT calculations to present, for example, lattice constants and elastic properties 
for pure Zr as well as some of the known hydride polymorphs [HI El [TOl HU [12]. Do- 
main et al. used ultrasoft pseudo potentials (UPP) and plane- waves single electron wave 
functions together with a GGA exchange-correlation (xc) functional PW91 to compute 
the structural parameters of the hydrides, the H-solution energies in Zr, and aZr elas- 
tic coefficients [Bj. In a second paper the same group presented DFT-calculations, using 
same methods as their first work, of stacking-fault energies and surface excess energies of 
solid solution H-a-Zr system [9]. Holliger et al. combined DFT-calculations, using UPP 
and plane-waves, similar to Domain's work, with cluster expansion (CE) in order to find 
the structures of several metastable H-Zr phases [TU]. Zhu's group presented in a recent 
study the structural as well as the elastic coefficient of the four ZrH x polymorphs also 
utilizing UPP and xc functional PW91. They, further, presented an analysis of the rela- 
tive phase-stability w.r.t temperature for the various hydrides using phonon-calculations 
evaluated at the gamma-point only [11]. Zhang et al. recently presented electronic struc- 
ture calculation on e-ZrH 2 using the projector- augmented wave (PAW) method together 
with xc-functional PBE[12j. They presented elastic constants and phonon dispersion 
and DOS curves. When comparing verious groups' results it is worth noticing that (a) 
Zr is notoriusly difficult to model and (b) all of the above cited works, except the last 
one, use a pseudo potential approach, which is generally considered less accurate than 
projector based and all-electron methods, such as PAW, FP-LAPW and FP-LMTO. In a 
recent paper Udagawa et al. even presented an analysis of plane defect in Zr-H systems 
by computing 7-surfaces and surface energy in these systems using DFT-methods[13j. 
There still exists some gaps and some cases where published elastic coefficients do not 
match from one group to another. Zhu et al. recently published a list of elastic coefficient 
of the four ZrH x hydrides but did not note any other calculations made of these [11]. In 




(a) aZr (b) 7 -ZrH (c) 5-ZrHi.g (d) e-ZrH 2 



Figure 1: The unit cells for aZr and the different hydrides used in the calculations 



the case of e-ZrH 2 , however, at least one other group have published a full set of elas- 
tic coefficients as well as phonon-based thermodynamical properties, and in the case of 
some of the elastic coefficients there are a clear discrepancy between the groups [12], 
even though both groups used similar DFT methods. A fourth hydride, called £-Zr 2 H, 
was recently detected by Zhao et al. [14] and will be part of future investigations. In 
the present study, however, only the most common hydrides (5-, 7-, and e-HZr) will be 
covered. 

In light of the importance of the elastic properties to understand stress driven embrit- 
tlement and the trouble with which to find experimental data on pure hydrides it is clear 
that DFT-calculations is a powerful tool. In order to bring clarity and fill in some blanks 
it is the goal of this work to first present the atomic and crystal structure of the various 
ZrH-polymorphs as well as the elastic coefficient. In the papers second part the phonon 
band structure and the free energies are presented, for the first time for some hydride 
phases. Accurate ground state electron structure calculations are essential in order for 
phonon-based thermodynamical properties to be accurate. From the phonon bands it 
is possible to compute the Gibb's free energy of formation as function of temperature 
and Helmholtz' free energy of formation as a function of pressure. These are necessary 
in order to properly deduce the phase diagram as well as needed in phase-field methods 
that are able to model the micro-evolution of the hydrides, for example as a function of 
temperature or load. 

Calculations 

DFT calculations are performed using the GPAW-code]15j. GPAW implements the 
PAW formalism using real-space grids to represent wave functions, pseudo electron 
densities and potentials. In these calculations the GGA type Perdew-Burke-Ernzerhof 
(PBE) exchange-correlation functional [TB], a real-space grid spacing of h=0.15 A and a 
Monkhorst-Pack k-point sampling[T7] with 6x6x6 k-points are used. 

The grid-spacing and number of k-points used in the calculations was chosen to ensure 
that the total energy was sufficiently converged with respect to parameter of interest. 
Increasing the number of k-points or decreasing the grid-spacing further was seen to 
change the total energy in the calculations with at most 0.05 eV. 



Equilibrium structures 

The equilibrium lattice parameters for a-Zr and the different hydrides (see figure [l] 
are determined by relaxation of ionic positions and calculation of the total energy of 
the compound in question for different values of the lattice parameters. A forth order 
polynomial is then fitted to the data and the minimum of this polynomial corresponds to 
the equilibrium lattice parameters. The 5-ZrH x is non-stoichiometric and experimentally 
found to have x=1.66. In order to model this phase we use the fee x=1.5 structure (cf. c 
in figure [T]). 7-ZrH was modeled using unitcell where hydrogen occupy the four tetragonal 
sites in the [110] plane. Another possibility exists, where four H occupy tetragonal sites 
in a diamond configuration. It is not included in the results of this work, however, since 
it has been found to have a total energy of formation about 0.3 eV higher than that of 
planar configuration. 



Elastic constants 

The elastic constants of the different zirconium hydrides are calculated by deforming the 
bravais lattice, R, of the crystal to a deformed state, R' according to: 

1 &xx 2^ X V 2&xz 

R = R -e 1 + e -e 
±k 2 x v yy 2 y z 

2&xz 2^y z ^ zz 

The change in energy of the crystal due to the deformation can be described by the 
following formula: 

i=l j=l 

This equation is valid for small strains and the coefficients is the elastic constants 
of the crystal and and e 3 - are components of the strain matrix that describes the 
deformation, where the indices i,j = 1,...,6 corresponds to the different strains, i,j = 
1, 6 = xx, yy, zz, yz, xz, yx. 

When the elastic constants for the material is known bulk modulus and shear modulus 
of the different hydrides can be calculated using standard relationship between B,G and 
E PS]. 



(1) 



Thermal properties 

The phonon band structure and thermal properties of the different hydrides are calcu- 
lated using the small displacement method as implemented in the phonopy[TP] code. In 
short, the free energy consist of terms involving the atomic energy, the enthalpy and the 
entropy. The first is simply the total energy provided by the zero-kelvin DFT calcula- 
tion. The enthalpy and the entropy both include temperature dependent terms which 
are functions of the phonon states (bulk vibrations). The phonons, in our calculations, 
was found by the finite displacement method whereby atoms were displaced a small 
distance. For the present study the distance was 0.1 A since it was found to be the min- 
imal displacement needed in order to accurately calculate the Hellman-Feynmann forces 



Table 1: The resulting equlibrium structures for aZr and 
the different zirconium hydrides together with experimen- 
tal data and the results from previous calculations 

present work experiment previous work 





a 


c 


a 


c 


a 


c 


a-Zr 


3.237 


5.157 


3.232 a 


5.147 a 


3.223 b 


5.175 b 








3.232 c 


5.148 c 


3.23 d 


5.18 d 


7-ZrH 


4.592 


4.998 


4.596 c 


4.969 c 


4.586 b 


4.9 b 












4.58 d 


5.04 d 


5-ZrHi. 5 


4.775 




4.771 a 




4.67 b 














4.79 d 




e-ZrH 2 


4.999 


4.433 


4.975 c 


4.451 c 


4.72 b 


4.21 b 








4.975 c 


4.47 e 


5.01 d 


4.44 d 



a Reference [20] , b Reference [TT] , c Reference [21] , 
d Reference [H] , e Reference [22] 



with amplitudes much greater than numerical error noise that would otherwise introduce 
chaotic behaviour in the results. A supercell of 2 x 2 x 2 unit cells matrix was used to 
calculate the phonon band structure for the different hydrides. For aZr a supercell with 
3x3x2 unit cells and a displacement of 0.1 A were used. 

Results 

The resulting equilibrium lattice parameters for the different hydrides and a-Zr is sum- 
marized in table [T] together with experimental data and the results from previous cal- 
culations. As seen the results from this study is in good agreement with experimental 
results and also with some of the previous calculations. 

0.14 I 




Figure 2: Total Energy relative to ground state energy vs. c/a for e-ZrH 2 

The results from the calculation of the elastic stiffness constants, bulk modulus, shear 
modulus and Young's modulus for the different hydrides are presented in table [2] together 



Table 2: Calculated elastic constants (C7y), Bulk modulus (B), shear mod- 
ulus (G), and Youngs modulus (E) for the different hydrides, everything in 
units of GPa. Comparison with experiments and previous calculations is 
presented if available 



phase 


method 


C n 


C33 


C\2 


Cl3 


C44 


Cm 


B 


G 


E 


aZi 


This work 


157 


158 


51 


62 


15 


44 


91 


29 


80 




PW91/UPP a 


142 


164 




64 


29 


39 


92 








PW91/UPP C 


160 


182 




66 


18 


51 


96 


39 


103 




PBE/UPP C 
Expt. b 


139 


163 


71 


66 


26 




93 








155 


173 


67 


65 


36 


44 


97 






7-ZrH 


This work 


131 


176 


123 


92 


64 


75 


117 


32 


88 




PW91/UPP C 


128 


187 


126 


70 


55 


66 


117 


44 


117 


5-ZrHi. 5 


This work 


108 


108 


139 


139 


51 


51 


128 


24 


69 




PW91/UPP C 


63 


65 


28 


44 


93 


101 


47 


63 


130 




Expt. f 














126 


50 


132 


e-ZrH 2 


This work 


156 


132 


144 


107 


40 


60 


125 


25 


73 




PW91/UPP C 
PBE/PAW d 


102 


108 


20 


11 


36 


24 


44 


37 


87 




166 


146 


141 


107 


31 


61 


130 


29 


80 



a Reference j8] , b Reference (23] , c Reference [IT] , d Reference [12] , 
e Reference [24J , f Reference [20J 



with results from previous calculations and experiments. The elastic constants of the 
present work agree nicely with the experimentally determined parameters, where they 
exist. It is also clear that for the e-hydride our results agree with those of Zhang et al. [12]. 
As can be seen in figure [2] we also obtain the ground-state geometry of c/a < 1 consistent 
with experiment. Previous calculations by Ackland|25j using Pseudo potential/plane- 
wave methods with PW91 xc-functional found that the lowest energy geometry was 
when c/a > 1. This further shows the difficulty of this system and supports the present 
results. The elastic properties of this hydride do not agree with those presented by 
Zhu et al. [11] however. It is not clear to us what is the cause of this discrepancy. 
Zr is known to be more difficult to calculate accurately than many other elements and 
since Zhu's group used pseudo-potentials to describe the core-electrons, and plane-waves 
for the valence electrons, it is possible that the complex Zr-H bonds are not described 
fully enough with their method. Our group and Zhang's, on the other hand, used PAW 
formalism based methods where the full all-electron wavefunction are fully described. 
Here the projectors are defined either on a grid in k-space or in real space in the case of 
Zhang's group (PAW) and ours (GPAW) respectively. 

The calculated phonon band structures for aZr and the different hydrides can be seen 
in figure [3j For gamma hydride the unit cell consists of 2 Zr and 2 H and for the delta 
hydride the phonon dispersion was calculated using a unit cell of 4 Zr and 6 H. The other 
hydrides and the hcp-Zr were calculated with primitive unit cells. The band structure 
for aZr has recently been calculated by Souvatzis and Eriksson [2S] arid their results are 
in good agreement with our calculations. The phonon band structure for e-ZrH 2 was 




Figure 3: The phonon band structure for aZr and the different hydrides, are frequencies 
are given in units of meV 



calculated by Zhang et al p2] and their results are in qualitative agreement with ours. 
This is is further support of the present results. 

In figure [3] the phonon density of states for aZr and the different hydrides can also be 
seen to the right of the phonon dispersion. In these figures we clearly see that there is a 
large gap in the phonon band structure for the different hydrides. We have low frequency 
vibrations from the heavy zirconium atoms and high energy vibrations from the hydrogen 
atoms. We also calculated the free energy of formation of the different hydrides, in table 
[3] the free energy of formation for the different hydrides at room temperature and 625 K 
( 350° C). From this table we can see the linear decrease of formation energy for every 
new H-Zr bond. Of course, this does not imply that e-ZrH 2 always have lowest formation 
energy, and thereby would always be preferred by nature relative the other hydrides; only 
that at hydrogen concentration that is stoichiometric relative to the particular hydride 
these values holds. It simply follows from the lower energy the system feel per added 
new Zr-H bond. To compare relative thermodynamical stability at a general hydride 
concentration one needs to formulate formation equation for each hydride in relation to 
a given hydride and then compute new formation energies. This will be done in future 
studies. 

Table 3: The Free energy of formation (in eV) for the three different hydrides at room 
temperature and 625 K ( 350° C) and at zero pressure. 





295 K 


625 K 


7-ZrH 


-0.6002 


-0.3528 


tf-ZrHi.5 


-0.8820 


-0.5240 


e-ZrH 2 


-1.2101 


-0.7214 



Conclusion 

The structural and thermodynamical properties of the three most common zirconium 
hydride polymorphs have been presented. Based on a novel all-electron dft-method 
using PAW formalism in evaluated using a real space grid sampling we have presented 
lattice parameters, elastic stiffness constants at zero pressure and temperature. In order 
to learn more about the finite temperature effect the phonon dispersion and DOS was 
calculated. Thermodynamical properties at two temperatures where finaly presented. 
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